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ABSTRAGT 

We present a technique to construct a spectropolarimetrically accurate magneto-hydrostatic model 
of a large-scale solar magnetic field concentration, mimicking a sunspot. Using the constructed model 
we perform a simulation of acoustic wave propagation, conversion and absorption in the solar interior 
and photosphere with the sunspot embedded into it. With the 6173A magnetically sensitive photo- 
spheric absorption line of neutral iron, we calculate observable quantities such as continuum intensities, 
Doppler velocities, as well as full Stokes vector for the simulation at various positions at the solar disk, 
and analyse the influence of non-locality of radiative transport in the solar photosphere on helioseismic 
measurements. Bisector shapes were used to perform multi-height observations. The differences in 
acoustic power at different heights within the line formation region at different positions at the solar 
disk were simulated and characterised. An increase in acoustic power in the simulated observations of 
the sunspot umbra away from the solar disk centre was confirmed as the slow magneto-acoustic wave. 
Subject headings: Sun: magnetic fields - Sun: oscillations - Sun: helioseismology - sunspots 


1. INTRODUCTION 

Techniques of local helioseismology are currently un¬ 
able to unambiguously determine sub-surface structure 
of the flows and sound speed perturbations in and around 
large-scale solar mag netic field concentrat i ons, such as 

a ots and pores (jShelvag et al.l l2QQ7at iGizon et al.l 
Rioradi et al.ll2QIQ[ ). Due to the complexity of mag¬ 
netohydrodynamic processes involved, our understand¬ 
ing of the behaviour of magnetoacoustic waves as they 
are absorbed, reflected and refracted by sunspots is far 
from complete. 

Recently, it was demonstrated that solar magnetic 
fields and the process of magneto-acoustic wave mode 
conversion associated with them lead to significant 
changes in the wave travel tirnes u s ed in helioseis¬ 
mic inversions (jMoradi fc Callvl l2Ql4 iHansen fc Callvl 
l2QI4f ). Eour processes associated with strong magnetism 
domi nate wave behaviour in sunspots (jCallv et al.l 
|2QI5[ ): fast/slow mode conversion at the Alfven/acoustic 
equipartition level Va = Cg, allowing acoustic (slow) 
waves to transmit into the upper atmosphere if the ‘at¬ 
tack angle’ a between wave vector and magnetic field is 
small b ut converting them to magnetic (fast ) waves oth¬ 
erwise (jCallvl [2QQ5 iSchnnker fc Callvl l2QQ6f ): the “ramp 
effect” that reduces the effective acoustic cutoff frequency 
ujc to ujc COS 0^ wh ere 0 is the magnet ic field inclination 
from the vertical (jBel fc LerovllI97^ : fast wave reflec¬ 
tion around the height where the Alfven speed matches 
the wave’s horizontal phase speed; and fast/Alfven 
mode conversion that typically occurs over several scale 
heights near the fast wave reflection level, generat¬ 
ing bo th upward and down ward propagating Alfven 
waves (jCallv fc HansenI l2QIIf ). East/slow conversion is 
found to produce large negative travel time shifts, while 
fast/Alfven conversion generates countervailing positive 
shifts provided the vertical plane containing the wave 
vector is nearly perpendicular to the vertical plane con - 
taining the magnetic field lines (jCallv fc Moradill2QI3f ). 


The ramp effect allows field-guided acoustic waves to 
enter the atmosphere in inclined field where frequency 
uj > uJc cos 0 while normally the acoustic cutoff would 
prevent their propagation {uj < cjc)- The complexity and 
sensitivity to magnetic field direction of wave motions 
above the equipartition level makes interpretation of ob¬ 
servations difficult but potentially rewarding. Moradi 
et. al. (2015, submitted) studied the effects of direc¬ 
tional time-distance helioseismology on the travel time 
measurements in the sunspot model. 

There are also possible significant discrepancies in 
travel time measurements originating from the effects 
of non-locality of radiative transport in the solar at¬ 
mosphere. Changes in spectral line forn iation heights 
due to magnetic field presence (see e.g. iShelvag et al.l 
1200^ . systematic centre-to-limb variations in absorp¬ 
tion line formation (jShelvag fc Przvbvlskil[2QT^ . as well 
as instrumental effects, such as stray light, and other 
processes involved in formation and measurement of ra¬ 
diation intensities and Doppler shifts result in our in¬ 
ability to unambiguously measure the travel time pertur- 
batioi^ and, the refore, infer solar sub-surface structure 

Rapid improvements in computational power already 
make it possible to perform forward modelling of mag¬ 
net ohy dr o dynamic wave propagation and mode con- 
version in “realistic ” solar magne t ic fie l d structures 
(IShelvag et al.l I2QQ9I: iMoradi et al.l l2QQ^ lEelipe et al.l 

2QIQI: iCameron et al.l I2QIII: iKhomenko fc Callvl 201 21: 
EelipelUQI2l: I Zharkov et al.l I2QI3I: lEelipe et al.l l2QI4f) . 
Spectral line synthesis codes and radiative diagnostics 
tools also allow computations of mock observables from 
the simulated plasma parameters, allowing for direct 
comparison between simulations and observations in 
computational helioseismology. 

Creating a sunspot that is both spectropolarimetri¬ 
cally accurate and magnetohydrostatically stable is in¬ 
herently difficult, as the sound speed and temperature 
can change significantly with small changes in the den- 
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sity and pressure strat i fcatio n. The sunspot model of 
iKhomenko fc ColladosI (j 2 QQ 8 f ) was created to allow em¬ 
pirical quiet and umbral solar models to be used in 
the near-surface layers in c ombination with a Schluter- 
Temesvary flux tube model (jSchluter fc Temesvarvlll958f) 
in the interior. However, the model created this way 
is still not convectively stable. Convective instability 
is fatal to linear MHD simulations, but these codes are 
less expensive than full non-linear simulations, and ideal 
for the long time series required in helioseismology; for 
the study of fast and slow magneto-acoustic waves; and 
for simulating fast-slow and fast-Alfven mode conver¬ 
sion in the photosphere and lower chromosphere. The 
effects of convective stabilisation on the eigenmodes of 
solar models for helios eismic simulations were studied by 
iSchunker et al.l (j 2 Qllf ). A technique for stabilising the 
atmosphere is discussed in Sec. O 

In this paper, we present a model of a magneto¬ 
hydrostatic and spectropolarimetrically accurate 
sunspot. Our model is based on t he sunspot-like 
model of IKhomenko fc ColladosI (j 2 QQ 8 f ). The model 
was adjusted to provide a more accurate replication 
of photospheric sunspot properties taken from semi- 
empirical models, while still maintaining a smooth 
transition of physical properties between the magnetic 
and non-magnetic regions required for stable numerical 
simulation. This technique makes it possible to obtain 
accurate photospheric absorption line formation heights 
as well as allowing the study of observational signatures 
of acoustic wave propagation in the simulated model 
at different positions on the solar disk. We perform a 
magnetohydro dynamic simulation of the propagation of 
a wave through this sunspot-like model and investigate 
the behaviour of acoustic waves in the simulated model 
using the synthesised radiation, as if it were observed. 
We also investigate effects of the centre-to-limb variation 
effects on Doppler velocity measurements and study 
the line bisector shapes to allow for a multi-height 
view in the line formation region, which can be used 
to observationally disentangle wave mode conversion 
process in the solar atmosphere. 

The structure of the paper is as follows. In Section 2 we 
describe the background model. In Section 3 we explain 
the magnetohydrodynamic simulation and the spectral 
synthesis methods used to provide artificial observables. 
Section 4 provides results and description of the radiative 
effects on acoustic wave measurements in observations. 
In Section 6 , we discuss our findings. 

2. MODEL 


where the gravity acceleration g, Brunt-Vaisala fre¬ 
quency N, and the adiabatic index 7 are functions of 
depth. The equations are solved using a fourth-order 
Runge-Kutta method on a one-dimensional grid. The 
equispaced grid covers the height range from —50 to 
+2.48 Mm and is resolved with the vertical step Az = 
0.02626 Mm. 

To enforce convective stability, the Brunt-Vaisala fre¬ 
quency dependence on z is modified by setting the neg¬ 
ative values in the convectively unstable solar interior 
to small positive ones. The free parameter a must be 
greater than zero to enforce convective stability and is 
increased so to match the p ressure in the model with the 
Stand ard Solar Model S (jChristensen-Dalsgaard et al.l 
pressure at a point in the interior z = —5 Mm. 
The quiet Sun model is created by smoothly joining the 
VALIIIC (jVernazza et al.lflOMl ) to Model S and integrat¬ 
ing from the top of the solar atmosphere downwards. 
Using the above, a convectively stable quasi-solar model 
can be created with negligible change to the photospheric 
line formation regions and the correct sound speed pro¬ 
file below the surface at the expense of a slightly reduced 
density and pressure in the deeper regions. 

Following _ the met hod described by 

IKhomenko fc ColladosI (j 2 QQ 8 f ). an axisymmetric sunspot 
model is created in cylindrical geometry using three 
parameters a, 77 and Bq which change the sunspot 
radius, magnetic field inclination and strength, re¬ 
spectively. A full description of the effects of these 
pa rameters on the magnet ic field configuration is given 
by IKhomenko fc ColladosI . The model is defined on a 
two-dimensional r-z plane discretised into a domain 
from —10 to 2 Mm in height, with a radius of 100 Mm 
and resolution of Az = 0.1 Mm and Ar = 0.2 Mm. 

Below — 1 Mm depth a Low-type magnetic flux tube 
(|Lowlll980[ ) is constructed using a n extension of the exact 
Schlu ter-Temesvary formulation (jSchluter V TemesvarvI 

For the near-surface layers z > —1 Mm and in t he at¬ 
mosp here a Pizzo-type magnetic flux tube is used (jPizzol 
Il986f ). The Pizzo method creates a pressure-distributed 
magnetic field structure through an extension of the Low 
formulation used above. The magnetohydrostatic equa¬ 
tion for a non-twisted, cylindrical structure can be sim¬ 
plifie d by introducing a magnetic vector potential (jLowl 
Il975f ) . This allows the probl em to be red uced to a single 
equation for a scalar u{r^z) (jPizz 311911 ) 

1 du d^u 2 dp{u, z) 

+ ( 3 ) 


To provide a convectively stable quiet Sun 

background model, th e _ method described by 

iParchevskv fc Kosoviche'^ (l2nnl is used. The Brunt- 
Vaisala frequency N‘^ = ^ ~ ^ must be positive for 

convective stability. Rearranging this equation in terms 
of dpjdz^ combining with the equation for hydrostatic 
stability (Equation [1]) and introducing a free parameter 
a gives. 


dp gp^ pN‘^ 

r\ ^ 1 

dz 7 p g 


where p{u, z) is the gas pressure along the field lines. 

The Pizzo method boundary conditions require both a 
quiet Sun (denoted with index q) and umbral (denoted 
with index um) pressure, density, temperature, and pres¬ 
sure scale height {h = ^) and temperature distributions 

as functions of depth. The quiet Sun model {pq^pq^hq) 
generated above was used for the outer boundary condi¬ 
tion. F or the inner boundary the Avrett semi-empirical 
model (jAvrettlll98l[ ) is used, which is then joined to the 
pressure and density profiles at the axis of the self-similar 
flux tube using log-linear interpolation. This is then con¬ 
vectively stabilised using Equations m and o as de¬ 
scribed for the quiet Sun above. The Wilson depression 
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can be prescribed by shifting the log(r5ooo) = 0 of the 
umbral model {pum^ Pum^ hum)- The pressure and scale 
height are then distributed throughout the domain using 
the following: 

p{u,z)=Pg{z) - {Pg{z) - Pum{z)) (^1 - >( 4 ) 

h{u,z) = hg{z) - {hg{z) - hum{z)) (l- ) (5) 

V u{nr,0)J 

The potential solution given by Equation m is used as 
an initial guess. The pressure distribution given by Equa¬ 
tions 0 and © is iterated together with Equation © 
using a Gauss-Seidel method. Thus, the complete force 
balance is calculated with a specified precision, giving a 
final distribution of the potential and pressure. 

The Pizzo and Low type flux tubes are then joined 
at z = — 1 Mm and recalculated using Equations iEj). 
The density and the radial and vertical components of 
the magnetic field vector 5^, Bz are calculated according 
to: 


_ p{r, z) 

9 iz)hir,z) 

(6) 

, . 1 du 

(7) 

Brir,z) = — — 
r oz 

, . 1 du 

(8) 

Bz{r,z) = - — . 

r or 


To extend this model below 2: = —10 Mm a vertical 
flux tube with a constant Bz and zero Br is used, and 
the pressure and density profiles are continued smoothly 
downwards. 

Einally, the EreeEOS equation of state (jirwini l2Q12f ) 
is applied to find the adiabatic index, temperature and 
sound speed at each grid cell in the model. The model is 
then converted to Cartesian geometry, giving the full set 
of physical parameters required for the MHD simulations 
and radiative transfer calculations. 

Using the procedure explained above the magnetic field 
structure pictured in Eigure [T] was constructed. The 
background image in the figure shows the modulus of 
magnetic field B. The field lines are nearly vertical in 
the “umbral” region (r < 10 Mm), and show inclination 
of about 60° in the “penumbral” region, r > 10 Mm, of 
the sunspot model. 

In the figure, the dashed line shows the log(r5ooo) = 
0 layer, while the dotted contours represent Cq/va = 
0.1, 1, and 10 levels. As is evident from the figure, 
in the umbral region at the axis of the sunspot, the 
log(r5ooo) = 0 layer is positioned higher than CsjvA = 1 
layer, suggesting formation of the continuum radiation 
in the magnetically-dominated sunspot atmosphere. 

The 6173A photospheric absorption line of neutral iron 
is used for observations of the full solar disk with the 
Helioseismic Magnetic Imager (HMI) onboard the So¬ 
lar Dynamic Observatory (SDO). Therefore, this line 
was chosen to carry out radiative diag nostics of the 
sunspot model using the SPINOR code (jSolankil [l987l: 
iShelvag et al.l l2QQ7b( ). Eor each one-dimensional col¬ 
umn of the model, continuum intensity and spectral line 
profile calcu lations are p erformed by solving the Unno- 
Rachovsky (jUnnol [T956f ) radiative transfer equation for 



Radius [Mm] 

Fig. 1. — Magnetic field structure of the sunspot model. The 
magnetic field strength is shown with the magnetic field lines over¬ 
plotted (solid). Also shown are the CsjvA — 1 (middle dotted), 0.1 
(upper dotted) and 10 (lower dotted) contours. The dashed line is 
the log(r5ooo) = 0 contour, representing the visible photosphere. 
Note that the aspect ratio is severely stretched. 

the Stokes vector / = [/, U, Q,^]. Off-disk centre ob¬ 
servations are simulated by inclining the numerical do¬ 
main and interpolating the density, temperature, mag¬ 
netic field and velocities onto the new line of sight {los). 
The slanting is performed around the z = 0 km height 
and in the direction of positive y (Eigure [2]). The ve¬ 
locity and magnetic field vectors are then projected into 
the new reference frame. The calculation uses 500 wave¬ 
length points with a = 0.002A to ensure the spec¬ 
tral line is highly resolved. The los velocity is given by 
Uos = cos6> -h Vx sinO. The magnetic field is recalcu¬ 
lated using a similar relation. 

Eigure [2] shows the continuum images of the sunspot 
model calculated for 0°, 30° and 60° angles between the 
surface and the los, which correspond to viewing cosine 
// = !., 0.866 and 0.5, respectively. We find that the 
model produces a realistic limb darkening dependence 
with a continuum value of 79% of the disk centre intensity 
at /i = 0.5. This is only slightly higher t han the 75% 
of the limb darkening curve determined bv lEoukal et al.l 

(I200I . 

The velocity response functions of the 6173A spectral 
line are shown in Eigure [3] for the quiet Sun, two penum¬ 
bral regions at ±10 Mm, and in the centre of the sunspot 
umbra for the chosen positions at the solar disk. These 
locations have been marked with crosses in Eigure O 
Since, for observations away from solar disk centre, two 
points at the same distance from the sunspot axis are 
not equivalent, the penumbral models have been chosen 
so that los of PI crosses the umbral region, while the los 
of P2 inclines further into the penumbra. The response 
functions were calculated by computing a perturbed pro¬ 
file with a small positive (directed towards the observer) 
los velocity perturbation and subtracting from it an un¬ 
perturbed profile for the same location. 

The top row of Eigure [3] shows the response functions 
of the four points in the model for the disk centre cal¬ 
culation. The perturbation is directed towards the ob¬ 
server, causing the line to be blue-shifted. The lobes of 
the response function tilt inwards towards the line core. 
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Fig. 2.— 6173A continuum intensity of the sunspot model at 
the observational angles (top to bottom) 0°, 30° and 60° to the 
vertical. The figures have been normalised to the quiet Sun value 
at 0° inclination. Inclination is performed towards an observer 
displaced in the negative y direction. The crosses show the two 
penumbral, and one umbral point used in Figure [3] 

marked as 0.0A in the plots, clearly demonstrating de¬ 
pendence of the sensitivity of the profile on height; the 
regions closer to the line-core (on the x-axis) are formed 
higher in the atmosphere. 

The top-right figure shows a fully Zeeman-split pro¬ 
file in a strong umbral magnetic field. Notably, while 
the line formation height range is narrower compared to 
the quiet Sun, the response function lobes are wider in 
wavelength, suggesting higher sensitivity of the line to 
velocity perturbations. 

The two penumbral points are identical in the so¬ 
lar disk centre simulation due to the symmetry of the 
sunspot. As the penumbral magnetic field is weaker, the 
line is not completely split. In the line core, the response 
function shows two smaller regions of sensitivity to the 
velocity perturbation. 

Figure [3] demonstrates that the line formation height 
range increases with the inclination angle. In the case of 
the quiet Sun (left column of the figure), it increases from 
^ 400 km at the disk centre to ^ 800 km at // = 0.5. As 
the line width does not change significantly, the wave¬ 
length range of the response function does not change 
with the inclination. 

For the cases of magnetised penumbral and umbral at¬ 
mospheres, the observed visible sunspot surface increases 
with the inclination angle. Between fi = 1 and 0.5, the 
logTsooo = 0 level for the penumbral points PI and P2 
is shifted downwards by ^ 100 km, and by ^ 400 km 
for the umbra. The line sensitivity height range also in¬ 
creases further away from the disk centre, similarly to 
the quiet Sun. 

Notably, the line profiles and the response function 


shapes for PI and P2 are very different. The far-side um¬ 
bra (second column of Figure O PI in Figure [2]) will have 
a formation range that extends into the highly magnetic 
umbra. This can be observed as an increasingly split 
profile as inclination increases in the 2nd and 3rd rows. 
The near-side penumbral pixel (third column, P2) will 
similarly form in a region of lower magnetic field. Due 
to the inclination of the magnetic field, PI will measure 
a higher magnetic field strength along the los while P2 
will incline against the direction of field line inclination. 

The angle of the two ridges is seen to be larger in the 
umbral distributions than the quiet sun. For a small 
wavelength range in the quiet sun up to 500 km of the 
atmosphere will be measured. Comparitively, in the 
umbral distribution a similar filter would only sample 
around 100 km. 

The large range of responses seen in the different 
penumbral and umbral positions will lead to larger un¬ 
certainty in the observation height of velocity measure¬ 
ments. The impact of Zeeman-split profiles on velocity 
measurements is only amplified at higher inclinations. 

3. NUMERICAL SIMULATIONS 

We perform a simulation of acoustic wave propagation 
in the simulated model with the SPARC code. The code 
was designed to solve the linearised ideal MHD equa¬ 
tions for wave propag ation in a stratified solar environ¬ 
ment (jHanasogell2Qlll ). The version of the code we em¬ 
ploy for the simulations uses Message Passing Interface 
(MPI) to parallelise the computation and reduce compu¬ 
tation time. It uses an implicit compact 6th order finite 
difference scheme applied to the horizontal and vertical 
derivatives. An explicit filter is used to prevent numeri¬ 
cal instabilities in the solution. The boundary conditions 
used in t he simulation include a Perfectly Matched Layer 
(PML) (jHanasoge et al.l I2Q1Q[ ) at the top and bottom 
boundaries, allowing for efficient absorption of the out¬ 
going waves. A 12.5 Mm ‘sponge’ type absorbing layer 
is used on the side boundaries, which a dds a linear fric - 
tion term to the governing equations (iColoninsI I2QQ4D. 
The c ode includes a Lorentz force limiter (jRnmpel et al.l 
l2QQ9f ). which is required due to the high Alfven speed 
in the solar atmosphere. The limiter, although unphys¬ 
ical, prevents reduction of the time step and excessive 
computational times. The Alfven speed limiter is set 
at va = 125 kms“^, which is sufficiently high to al¬ 
low fast MHD waves of interest, which have horizontal 
phase speed less than this, to propagate and refract cor¬ 
rectly while still allowing us to use a manageable time 
step. Our cap is large enough for this to not be an oner¬ 
ous constraint. The implications of the limiter on helio- 
seismic travel time shi fts have been studied in detail by 
iMoradi fc Callvl (j2Q14f ). 

The numerical grid has horizontal extent of rix = riy = 
256 grid points, with a physical size of 140 Mm, giving 
a resolution of Ax = Ay = 0.55 Mm in the horizontal 
directions. To deal with the large variation in physical 
parameters over the domain from the convective zone 
to chromosphere, the code uses a vertical grid spacing 
based on the sound speed. The grid has Uz = 300 points 
between 1.5 and —25 Mm and is distributed such that 
the acoustic travel times between each cell are the same 
for the quiet Sun, Az oc l/cg. This gives a resolution of 
around 50 km near the photosphere, and around 1 Mm 
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Fig. 3. — los velocity response functions of Stokes / profile of 6173 A Fel line computed for the models of quiet Sun (first column), near-side 
and far-side penumbrae (second and third columns, respectively), and umbra (fourth column) for jj, = 1., 0.866 and 0.5 positions at the 
solar disk. The y-axis represents the height along the los at which the perturbation is placed, where 0 Mm represents the logrsooo = 0 
layer for the quiet Sun photosphere. The SOOOA optical depth axis is also shown, with a dashed line showing the observed depression of 
the photosphere. The corresponding Stokes-/ (top right) and Stokes-V (top-left) profile shapes were over-plotted in white in each panel 
over the wavelength range shown in the figure. 


in the lower solar interior. This means we do not resolve 
slow waves in the large /3 regime, but these are effectively 
decoupled from the system anyway, so their neglect is not 
important. The following a c oustic source, similar to that 
described by iShelvag et all (j2QQ9[) , was used: 


. . .27rt. —(t —Ti)^ — (f —ro)^ . . 

Vz = Ao sm( —) exp-^-2- (^) 


err 


X exp 


-{z - Zof 


where To = 300 s, Ti = 600s, at = 100 s, axy = 1 Mm, 
az = 0.25 Mm. The position of the pulse is ro{x,y) = 
(45, 70) Mm, zq = —0.65 Mm. 

The SPARC code solves the MHD equations for the 
perturbations around the MHS background model. A 
master-slave Open-MPI code has been written to take 
these perturbations, combine them with the background 
model and incline them as required. The SPINOR rou¬ 
tines are then applied to each pixel to generate the full 
Stokes vector for each pixel. Using the generated Stokes- 
/ profiles, the los centre-of-gravity Doppler velocity is 
calculated by computing the position of the centre of 
gravity of the line profile and determining its shift from 
the unperturbed counterpart, computed for the back¬ 
ground model, according to: 


= /(U-ZMA <“> 

To calculate bisector Doppler velocities from the spec¬ 
tral line the relative intensity Irei was determined by 
normalising the measured Stokes I between 0 and 1. Bi¬ 
sectors of the spectral line were calculated at 100 evenly 
spaced values between 0.05 — 0.95 of AeZ- The bisectors 
were calculated for the background model and for each 
output snapshot. A Doppler velocity was then deter¬ 
mined for each snapshot using the shift from the unper¬ 
turbed background value, according to: 

e^bsr — (Aq ^bsr^~{ • 

Ao 


4. RESULTS 

Using the model described in Section 2 and methodol¬ 
ogy given in Section 3, a 2.5 hour simulated observation 
of wave propagation through the sunspot model was per¬ 
formed. The top panel in Figure H] shows a time-distance 
plot of the centre-of-gravity los Doppler velocity mea¬ 
sured using Equation m- The first three wave bounces 
can be easily seen. A shift in the wave arrival time can 
be observed as a flattening of the wavefront as it passes 
through the sunspot umbra at ^ = 0 Mm. The middle 
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Fig. 4. — Response of the model to the acoustic source. Top 
panel - simulated Doppler velocities at 0° inclination, measured at 
X = 0. Middle panel - simulated vertical component of velocity at 
a geometric height z = 0 Mm, measured at tr = 0. The differences 
between these two can be seen to be small. Bottom pannel - ^/pVy 
through spot centre (0,0) Mm showing the propagation of slow 
modes down through the box. Two dashed horizontal lines in the 
top plots mark the position of the sunspot umbra. 


panel of Figure |4] shows a time-distance plot of the ver¬ 
tical component of velocity at the z = 0 Mm level of the 
simulation domain. A comparison between the top two 
panels shows that the vertical velocity in the domain and 
the los Doppler velocity are visually identical. Some re¬ 
flection can be seen from the top and bottom PMLs, and 
from the side boundaries. 

The bottom panel of Figure |4] shows the horizontal 
component of velocity, scaled by y/p to provide a view 
of the slow magnetoacoustic wave in the strong magnetic 
field. A slice is taken through the centre of the simu¬ 
lated sunspot (x = 0,p = 0). The fast wave can be seen 
to propagate through the sunspot in the lower interior 
where plasma-yd is high. At around z = —0.400 Mm in 
the umbra, the incoming fast wave hits Cs/va = 1 layer 
(See Figure [T]) and undergoes partial transmission as a 
slow mode (effectively acoustic in Cg < va)- The slow 
magnetoacoustic wave (now magnetic in Cg > va) can 
be seen to propagate back down into the sunspot as a 
flattening banding in the time distance plot. The wave 
amplitude in the atmosphere is low due to scaling by 
the very low densities, however, it still can be seen to 
continue to travel upwards above the photosphere and 
escapes through the absorbing upper boundary. 

Figure [5] shows a power spectrum plotted with az- 
imuthally averaged wavenumber and frequency. The 
spectrum has been calculated from the 2.5 hour time- 
series of centre-of-gravity velocities calculated from the 
synthesised line profile. The power ridges are well re¬ 
solved, although there are gaps in power at 4.5 mHz and 
shifts see n for high 1. Similar gaps are foun d in the simu¬ 
lations of iParchevskv fc KosovichevI (j2QQ7l ) with high top 
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Fig. 5.— A v — I power spectrum for the sunspot box calculated 
from synthetic Doppler velocities. 


boundary (1.75 Mm; their Fig. 6c), and attributed to 
trapping of acoustic modes. We do not understand how 
acoustic trapping explains this phenomenon. The gaps 
are also present in quiet sun simulations (no magnetic 
field), but are largely removed when the top boundary 
is lowered to 500 km above the solar surface (their Fig. 
6b). This suggests that there is some numerical dissi¬ 
pation mechanism operating in our model chromosphere 
{z > 500 km) that we are yet to identify. 

Using the Fel 6173 A spectral line data, velocity shifts 
were calculated for each bisector as the differences be¬ 
tween the perturbed and background values (Equation 
(HID). From the 2.5 hour time series of these velocities 
acoustic power maps of wave propagation in the sunspot 
were created. The power maps are calculated for 100 
bisector depth positions by taking the fast Fourier trans¬ 
form of the velocity time series for each of the (x, y) pixels 
in the model. The acoustic power is binned into differ¬ 
ent frequency bins by applying a Gaussian filter with a 
FWHM of 0.5 mHz around the chosen central frequency. 
For each frequency band and inclination the acoustic 
power is normalised by its average in the simulation do¬ 
main. 

By taking the Doppler velocity measurements at mul¬ 
tiple bisector depths, it was possible to disentangle the 
dependence of the wave behaviour on the height within 
the line formation region. As Figure [3] shows, bisectors 
taken higher in the line profile are formed deeper and 
closer to the continuum formation height, at a physical 
height of around 150 km, while bisectors taken deeper 
are closer to the absorption line core and formed higher 
at a physical height of around 500 km. As was already 
noted, these formation heights vary significantly in the 
penumbra and umbra of the sunspot, with an offset of 
350 km due to the prescribed Wilson depression. 

First we aim to investigate the horizontal distribu¬ 
tion of acoustic power throughout the simulation do¬ 
main. The region of the acoustic source at x = 45 Mm, 
^ = 70 Mm has been masked in the power maps. Fig¬ 
ures and [3 show the acoustic power calculated from 
the Doppler shifts measured at the bisector positions of 
0.3 Irei and 0.9 AeZ, respectively. The acoustic power 
was binned into frequency 1 mHz bands centred at 3, 4, 
5, 6 and 7 mHz (left to right in the figures), and the disk 
positions of —60°, 0° and 60° were used (top to bottom). 
los velocity measurements off the disk centre are affected 
by a larger line formation region and the presence of hor- 
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izontal velocity components in the los velocities. Imme¬ 
diately obvious in the figure is a series of concentric rings 
of power travelling out from the source. These rings can 
be thought of as the spatial analogue of the ridges seen 
in Figure [5] and occur at discrete values because a single 
point-like source was used in the simulation. 

The differences between the acoustic power maps in 
Figure [6] and [3 represent changes in the wave behaviour 
over the formation height of the spectral line. This for¬ 
mation height can span almost a megameter close to the 
solar limb (see Figure [3j), and the differences are substan¬ 
tially more pronounced in the magnetic regions. 

In the acoustic power maps (Figure [6] and [7|) a solid line 
representing the sunspot umbra is over-plotted. Outside 
the sunspot umbra the sunspot shadow is observed at 
all frequencies. In the sunspot shadow, the power from 
the pulse has been absorbed or reflected by the magnetic 
field structure. This is most obvious at lower frequencies 
(left two columns). As frequency increases, the ridges 
of increased acoustic power can still be seen behind the 
sunspot. Interestingly, the magnetic field perturbs the 
concentric rings seen at 7 mHz, as the fast wave prop¬ 
agation speed and turning height change. Behind the 
sunspot in the 7 mHz band (right column), a small region 
of increased power can be seen between the two outer¬ 
most rings of the power ridges (around ^ = 30 Mm). As 
this is seen in high frequencies and as a ring around the 
sunspot, rather than around the source, it appears to be 
a far-side acoustic halo. Acoustic halos are seen around 
active regions as an excess of acoustic power compared 
to the quiet sun. Q 

Comparison of the weakly-magnetic regions of each col¬ 
umn in Figure [6] shows little variation with inclination, 
regardless of frequency. In these regions the propagat¬ 
ing fast wave dominates the simulated observations as 
there is little magnetic field to alter the fast-wave turn¬ 
ing point or to allow for mode-conversion process to take 
place. Inside the umbral region, the acoustic power maps 
for the disk centre case show very little power in 3 and 4 
mHz, and the power in the vertically-aligned oscillations 
is seen to be almost completely absorbed. As the incli¬ 
nation increases, the 3 and 4 mHz power in the umbra 
remains low, while the 5, 6 and 7 mHz power bands show 
a significant power enhancement. 

From the response function (first coiumn in Figure [3j) 
we expect there to be iarger differences in power between 
the fine core (Fig. [6]) and fine wings (Fig. [7j) as we ob¬ 
serve further away from the soiar disk centre. There is 
iittie difference found in the disk centre cases (top row) 
between the two figures. However, at 60° inciination sig¬ 
nificant differences can be seen between the power maps 
at the fine core and fine wings. Particuiarly, (1) the um- 
brai power increase is only seen at the line core (Fig. [6j); 
(2) the ring-like structure (marked by the dashed circle 
in the left panel of Fig. El) found at around ^ = — 20 Mm 
in the bottom left two panels is somewhat wider at the 
line core than in the line wings (Fig. (3). This struc¬ 
ture is most apparent in the 3 — 4 mHz frequency bands 

^ A more comprehensive look at physics of acoustic halos, based 
on three -dimensional sim ulations of this sunspot model can be 
found in IRiis et ^ (I2015l f , which expands on the previous works 
of lHa.na.sogfil 1120081^ : IKhomenko Colladosl l|200®. They are at¬ 
tributed to fast MHD waves reflected in active region atmospheres 
by the steep Alfvn speed gradients there. 


(first two columns), and the power in it decreases with 
increasing frequency. While the inclination of the mag¬ 
netic field at the surface at the radius of 20 Mm is 60° de¬ 
grees from vertical the magnetic field strength is low, and 
the equipartition layer Cs = va is located above the line 
formation region. Therefore, the ring is of acoustic na¬ 
ture and cannot be related to the slow magneto-acoustic 
mode, as it is found at the source side in both —60° and 
60° inclinations corresponding to the los direction which 
is either parallel or highly inclined to the magnetic field. 

As demonstrated, the umbral and penumbral acoustic 
power structures are mostly seen near the line core (Fig¬ 
ure [6|). In the los velocities measured from bisectors in 
the line wings (Figure (3) only a faint structure can be 
seen at high inclinations, again more obvious in the 7 
mHz power band (bottom right). 

To better understand the three-dimensional structure 
of the umbral power increase, in accordance with the re¬ 
sponse functions shown in Figure [3l a multi-height obser¬ 
vation is made by computing the Doppler velocities from 
bisector shifts measured at different line depths within 
the line formation region. In Figure [HI a Doppler veloc¬ 
ity map for a slice marked by a dashed line in Figure [3 is 
plotted for each bisector depth in the range 0.3 —0.9 AeZ- 

For the 0° and 30° inclinations (top and middle rows of 
Figure IH)), the changes in the structure and magnitude of 
acoustic power over the line formation range are limited 
to an increase in power in the high frequency bands for 
observations made closer to the line cor e. This matches 
the o bservations of acoustic halos by iRaiaguru et ah! 
(HoH), where the acoustic power was weaker using filter- 
grams close to the line-wings. At 60° inclination (bottom 
row, left two columns) the faint x = ±20 Mm radius low 
frequency ring can be made out, increasing in radius at 
larger heights. 

The vertical extent of the umbral power structure, seen 
at high inclinations at 6-7 mHz (bottom right two pan¬ 
els in Figure IBHH)) shows a significant increase in power 
higher in the formation region. The formation of this 
acoustic power structure seems to start mid-way up the 
line formation region, suggesting a highly localised phe¬ 
nomenon. The inclination of the field lines at the centre 
of this power increase {y = 4.11 Mm) is approximately 
20° from the vertical. Taking into account the observa¬ 
tion angle of ±60°, the field line is almost perpendicular 
to the line-of-sight. Comparison with Fig. (3 shows the 
continuum formation height of the spectral line and the 
Cs/va = 1 layer cross at z = 0 Mm and 4 Mm radius 
from the sunspot centre, allowing for direct observation 
of conversion of fast (parallel to the magnetic field line, 
perpendicular to the los) waves to slow (perpendicular 
to the magnetic field line, and parallel to the los) waves. 
The increased power in this region corresponds to the 
slow magneto-acoustic wave in the region where the mag¬ 
netic field is close to perpendi cular to the los. This result 
confirms previous findings bv [Zharkov et al.l (j2Q13[ ) that 
the observed acoustic power increase in the sunspots is a 
signature of slow magneto-acoustic waves. 

5. DISCUSSION AND CONCLUSION 

In this paper we have: (1) described a modification 
of the Khomenko sunspot model we developed to pro¬ 
vide a more accurate line formation region, allowing for 
accurate spectral synthesis to be performed; (2) anal- 
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Fig. 6.— Acoustic power map calculated from the shifts in the bisector of the Fel 6173A line at a bisector height of 0.3 I^el- The columns, 
from left to right, show power in the 3, 4, 5, 6 and 7 mHz bands. The rows, from top to bottom, show measurements made at —60°, 0° 
and 60° inclination from the vertical, where the field of view has been inclined in the y-direction. The power at each inclination angle 
has been normalised by its average in each frequency band. This represents a velocity measurement made near the line-core, showing the 
behaviour higher in the atmosphere. The sunspot umbra has been marked with a solid circle, while the dashed line in the bottom right 
panel represents the slice taken in Figure [8] The low frequency ring has been marked in the top left panel. 
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Fig. 7.— Acoustic power map calculated from the shifts in the bisector of the Fei 6173A line at a bisector height of 0.9 Irel- The layout 
of the columns and rows is as seen in Figure [6] This figure represents measurements made near the line wings. 
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Fig. 8.— Bisector power map at x = 0 Mm. The layout of the columns is as in Figure[6] The rows, from top to bottom, show measurements 
made at 0°, 30° and 60° inclination from the vertical. The Y axis in this figure represents the line depth where the bisector wavelength 
and Doppler velocity are measured and covers a large part of the line formation region of 400 — 800 km in length. The formation region 
will depend on the magnetic field strength and inclination (Figure [3]). 
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ysed response functions in our model to understand 
our synthesised centre-to-limb observations of the Fel 
6173A spectral line in a model sunspot; (3) have investi¬ 
gated wave propagation through the model using a three- 
dimensional simulation of magnetoacoustic wave propa¬ 
gation; (4) computed the spectral line profiles and pro¬ 
vided a time series of the simulated observations of the 
sunspot model using the simulation data; (5) used spec¬ 
tral line bisector measurements to perform multi-height 
simulated observations over the line formation region; 
(6) studied maps of the acoustic power in the sunspot 
to better understand the absorption line response to the 
oscillations in sunspots and the effects of non-locality of 
radiative transport on helioseismic measurements. 

The sunspot is found to absorb or scatter the incom¬ 
ing acoustic wave energy in all but the 6 and 7 mHz 
frequency band. There are signatures of a slight power 
enhancement seen around the sunspot, similar to those 
seen in observations of acostic halos. 

Small ridges of acoustic power ca n be seen in the 7 mHz 
band in the shadow of the sunspot. I Zharkov et al.l (j2Q13f ) 
showed in a simple 2D simulation that there was a sig¬ 
nificant power enhancement as a far-side acoustic halo. 
We can see slight power enhancements in the range of 
25-40 Mm from the sunspot. Comparing the different 
rows in Figures [6] and [3 shows that outside the sunspot 
umbra, the horizontal and vertical velocities behave sim¬ 
ilarly, with only minor differences in the power maps. 

The appearance of a high-frequency anomalous acous¬ 
tic power excess in the sunspot centre - the umbral 
“belly button” - can be seen predominantly in the case 
of 60° inclination with faint traces at lower inclina¬ 
tions. This geometry suggests that it is driven by slow 
magneto-acoustic waves, as it is seen when the hori¬ 
zontal velocity component dominates the los velocity 
(bottom rows Figure [6] and (3). It can be seen as a 
crescent-like structure, which is most dominant towards 
the line core (higher in the atmosphere. Figure [6]) and 
very faint in the spectral line wings (lower in the at¬ 
mosphere, Figure [3 • Umbral p ower enhancem e nts ar e 
seen in the space- based (HML IZharkov et all (j2013l )) 
and ground based (jBalthasar et al.l Il998f ) observations 
of sunspots. Notably, no power exce ss was observed in 
the G-band (|Nagashima et al.l I2QQ7I ) . The appearance 
of this power increase in a simulated sunspot suggests a 
magneto-acoustic pheno menon, rather than photon noise 
(jPonea fc Lindsewl2Q15f ). There are many differences in 
both sunspot properties and the radiative effects on HMI 
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measurements that could explain the lack of the umbral 
power increase in all acoustic observations of sunspots. 
These include changes in the Wilson depression, the wide 
range of the velocity-response function in a magnetic 
structure (Figure [3|), low resolution and high-noise mea¬ 
surements off the disk-centre or issues with using discrete 
filters on highly split profiles. 

Current measurements of acoustic travel times in com¬ 
putational helioseismology are largely performed using 
measure ments at the geometric heights in the simulation 
domain (|Moradi fc Callvll2Q13F or on a surface roughly 
representing the continuum formation height determined 
by the log(r5nnn) = 1 layer in the siuinlatio n domain 
(|KhomenkQ fc CalK^ 12012t iMoradi et al.l I2Q15D . Despite 
the fact that the physical velocity in the simulation 
matches reasonably well to the los velocities calculated 
from the simulated spectral lines, this method misses a 
lot of information that can otherwise be gained from the 
range of formation of the spectral lines. As we show, 
this range also changes substantially if the simulation 
is performed for positions at the solar disk away from 
the centre. Using the model we described, artificial ob- 
servab les mimicking the HMI and MDI pipel ines can be 
made (jScherrer et al.|[2QT^ iFleck et al.ll2Qll[) . as well as 
comparisons to ground based observations. 

The multi-height D oppler measurements made by 
iNagashima et al.r(|2Q14D . using the HMI filter-grams pro¬ 
vide a similar approach to multi-height measurements as 
the bisectors used in this study. Rather than velocity 
measurements made using shifts in Stokes-/, HMI uses 
measurements of both Stokes / + U and Stokes I — V. 
As the next step, it will be important to fully simulate 
the HMI data pipeline response to a variable magnetic 
atmosphere before a direct comparison to the HMI mea¬ 
surements can be made. 
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